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Abstract 

In this review we discuss Cosmological N-Body codes with a special emphasis on Particle 
Mesh codes. We present the mathematical model for each component of N-Body codes. We 
compare alternative methods for computing each quantity by calculating errors for each of the 
components. We suggest an optimum set of components that can be combined reduce overall 
errors in N-Body codes. 

1 Introduction 

Observations suggest that the present universe is populated by very large structures like galaxies, 
clusters of galaxies etc. Current models for formation of these structures are based on the assump- 
tion that gravitational amplification of small perturbations leads to the formation of large scale struc- 
tures. In absence of analytical methods for computing quantities of interest, numerical simulations 
are the only tool available for study of clustering in the non-linear regime. The last two decades have 
seen a rapid development of techniques and computing power for cosmological simulations and the 
results of these simulations have provided valuable insight into the study of structure formation. In 
this paper we will describe some aspects of the cosmological N-body simulation, based on the code 
developed and used by the authors. We will stress the key physical and numerical ideas and detail 
some useful tests of the N-Body code. After an initial introduction to cosmological N-Body codes 
the discussion of N-Body codes will focus on Particle Mesh [PM] codes. [For a comprehensive 
discussion of numerical simulations using particles, see Hockney and Eastwood (1988).] 

An N-body code consists of two basic modules: one part computes the force field for a given 
configuration of particles and the other moves the particles in this force field. These two are called 
at each step to ensure that the force field and the particle trajectories evolve in a self consistent 
manner Apart from these we also need to setup the initial conditions and write the output. Thus the 
basic plan of an N-Body code follows the flow chart shown in figure 1. 

Structure of these modules depends on the specific application at hand. Here we outline some 
features that are particular to cosmological N-Body codes. To make quantitative statements about 
the required parameters of an N-Body simulation we shall assume that the particles populate a region 
of volume V. We shall also assume that we are using Np particles to describe the density field. The 
following physical requirements have to be taken into account while writing cosmological N-Body 
codes. 
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Figure 1: Flow chart for an N-Body Code. 
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The universe is filled with matter over scales that are much larger than the scales of interest 
in an N-Body simulation. Density averaged over larger and larger scales in the universe 
approaches a constant value. Thus the simulation volume V can not be assumed to exist in 
isolation and we must fill the region outside this volume by some method. Periodic boundary 
conditions are the only viable solution to this problem. In absence of periodic boundary 
conditions most of the matter in the simulation volume has a tendency to collapse towards 
the centre of the box, giving rise to a spurious, large rate of growth for perturbations. [A 
compromise solution, called quasiperiodic boundary conditions, has been tried to solve this 
problem for tree codes. In this scheme periodic boundary conditions are used to restructure 
the simulation volume to bring a particle at the centre of the box while computing force at 
its position. This is done for each particle so that the r e is n o strong tendency to collapse 
towards the centre of the box. ( Bouchet and Hemquisl ( 198^)] The most natural geometry 



for a volume with periodic boundary conditions is a cube. 

• We would like the evolution of perturbations to be independent of the specific boundary con- 
ditions. Thus the simulation volume should not be dominated by any one object. [If the 
simulation volume is dominated by one object then the tidal force due to its own periodic 
copies will influence its later evolution.] 

• We require the average density in the box to equal the average density of the universe. Thus 
perturbations averaged at the scale of the box must be ignorable at all times for the model of 
interest, i.e., a{R = V^^^) <^ 1. [Here cr{R) is the rms dispersion in the mass averaged at 
scale R.] For example, in case of the standard CDM model this would require the box to be 
at least 100/i^^Mpc [at z = 0] in extent. 

• We would like to probe scales that are sufficiently non-linear in order to justify the use of 
N-Body simulations. If we wish to compare the results with the real universe then we should 
be able to study the formation of structures with a large range in scales. In other words, the 
mass of individual particles in an ideal simulation must be less than the mass of the smallest 
structure of interest. Therefore the smallest number of particles required to cover the relevant 
range of scales to study galaxy clustering, say, is 

^ > — ~ 11 ~ 5 X 10^ (1) 

Mgalaxy 10^^^ 

We need a very large number of particles to represent the density field over this range of 
scales. Therefore most numerical techniques used in cosmological N-Body codes are oriented 
towards reducing the number of operations and stored quantities per particle. 

• We approximate collection of a very large number of particles in the universe by one particle 
in an N-Body simulation. Therefore the particles in an N-Body simulation must interact in a 
purely coUisionless manner. 

The periodic boundary conditions and a very large number of particles are the key considerations in 
developing the detailed algorithm for an N-Body code. Of the two main components in an N-Body 
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code, integration of the equation of motion is a process of order 0{Np). The calculation of force, 
if done by summing the force over all pairs, is a process of order 0{Np). Therefore the calculation 
of force is likely to be more time consuming than evolving the trajectories in N-Body codes with 
a large number of particles. It can be shown that, for particle numbers greater than a few hundred, 
direct computation of force takes an excessively long time even on the fa s test computers. [The ve ry 



high speed special purpose computers are an exception to this. (iBrieul jSummers and Qstrikeg) )1 
Three schemes have been evolved to address this problem and replace direct summation over pairs 
by methods that are less time consuming. These are 

• Particle Mesh (PM) : Poisson equation is solved in the Fourier domain and the potential/force 
is computed on a fixed grid. The force is then interpolated to particle positions for moving the 
particles. Density field, the source for gravitational potential, is also computed on the same 
mesh/grid from particle positions by using a suitable interpolating function. The "smooth- 
ing" of particles limits the resolution of such simulations but ensures coUisionless evolution. 
[See Bouchet, Adam and Pellat (1984) and Bouchet and Kandrup, (1985) for an excellent 
discussion of PM codes.] 

• Particle-Particle Particle Mesh (P^M) : This scheme (Efstathiou et al, 1985) improves the 
resolution of PM method by adding a correction to the mesh force for pairs of particles with 
separation of the order of, and smaller than, the grid length. The number of operations re- 
quired for this correction is proportional to Npn{R) where n{R) is the average number of 
neighbouring particles within a distance R. This can be written as Npn{l-\-^{R)), where n is 
the average number density and ^ is the averaged correlation function. It is obvious that such 
corrections can become time consuming in highly clustered situations [when ^{R) ^ 1]. 

• Tree : In this scheme, information about the density field is set up in form of a hierarchical 
tree. Each level of the tree specifies position of the centre of mass and the total mass for a 
region in the simulation volume. Force from particles in distant regions in the simulation box 
is approximated by the force from the centre of mass for particles in that region. This leads 
to a reduction in the number of operations required for calculating force. [Barnes and Hut 
(1986); Bouchet and Hernquist (1988)] To estimate the number of operations required for 
setting up the tree structure and evaluate the force, let us assume that the simulation volume 
is a cube and is subdivided into eight equal parts at each level. This subdivision of cells is 
continued till we have at most one particle per cell at the smallest level. Each particle is 
parsed at most once at each level, therefore the upper bound on the total number of operations 
is proportional to Np\n{L},ox/lm.in) where Imin is the smallest inter-particle separation for 
the given distribution of particles. We have, 

1/3 _ A-1/3^-1/3 _ ^-1/3 r r-1/3 ^2) 
"max ~ '^max "' ~ -^^p ^box'^max v^/ 



where n is the average number density, Smax is the maximum density contrast and rimax is 
the highest number density in the given distribution of particles. This implies that the upper 
bound on number of operations is proportional to 0{Nphi NpS^nax) ■ Incorporating periodic 
boundary conditions is a very nontrivial problem in the context of tree codes. [See Hernquist, 
Bouchet and Suto (1991) for a discussion of periodic boundary conditions for tree codes.] 
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Amongst these methods the PM scheme has two advantages over the other methods. Firstly 
it is the only one of the three methods outlined above that ensures a coUisionless evolution. [See 
Melott et al. (1997); Suisalu and Saar (1996).] Secondly it has the simplest algorithm and it is also 
the fastest of the three methods discussed above. In the remaining discussion we shall focus only 
on PM codes. However, apart from computation of force, all other components are the same in all 
these codes [with some minor variations] and most of the conclusions about relative merits of the 
available alternatives are applicable for all three types of codes. 



2 Moving the Particles 

In this section we will first present the system of equations that describe the evolution of trajectories 
of particles. This will be followed by a description of the methods used for numerical integration of 
equation of motion. 



2.1 Equation of IMotion 

It can be shown that the evolution of perturbations in a non-relativistic medium in an expanding 
background can be studied in the Newtonian limit at scales that are much smaller than the Hubble 
radius du = cHq^. The equations for a set of particles interacting only through the gravitational 
force can be written as 

a. 1 ^ 

X + 2-x = TjVx^ 

a 

Vlip = A7rGa^g5 = ^H^^o-- (3) 

Here the last equality follows from the Friedmann equations that describe evolution of the scale 
factor for the universe. The variables used here are : comoving co-ordinates x, time t, scale factor 
a, gravitational potential due to perturbations if, density g and density contrast 6. Cosmological pa- 
rameters used in this equation are : the Hubble's constant Hq and the density parameter contributed 
by non-relativistic matter Oq- 

N-Body simulations integrate this equation numerically for each particle. Numerical integra- 
tion can be simplified by modifying the equation of motion so that the velocity dependent term is 
removed from the equa tion o f motion. This can be achieved by using a different time parameter 9 
Jicivpin and ShandarinI (Il983h ). This is defined as 

de = Ho^ (4) 



a2 



In terms of which we have 

d^x _ 3 
de^ ~~2 



y\ = - (5) 

a 

In this form, all the cosmological factors can be clubbed in the source term in the equation of motion, 
making numerical implementation much simpler. 
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2.2 Numerical Integration 



In any N-Body code, a great deal of computer time is devoted to integration of the equation of 
motion. Tlie basic idea of time integration is simple : the equation of motion expresses the second 
derivative of position in terms of position, velocity and time. Time integration translates this into 
the subsequent changes in position and velocity. In the following discussion we will develop the 
idea of numerical integration of trajectories. 

Writing the position and velocity at time t + eina Taylor series about the position and velocity 
at time t, we have 



Here x is the position, v is the velocity, a is the acceleration and j is the rate of change of accelera- 
tion, and the subscript is used to identify the particle. We can use the equation of motion to express 
acceleration and the higher derivatives of position in terms of positions and velocities. Therefore, 
in principle, we can find the new position and velocity with arbitrary accuracy. However, such a 
scheme is not very practical from the computational point of view as it requires an infinite series 
to be summed. A more practical method consists of truncating the series at some point and using 
sufficiently small time steps to ensure the required level of accuracy. To illustrate this point, let us 
consider terms up to first order in the above series, we get 



This is the Euler's method of solving differential equations and the error in the solution is of order 
e^. Thus choosing a smaller time step leads to smaller error in each step. However, the number of 
steps required to integrate over a given interval in time is inversely proportional to the time step so 
that the overall error changes only linearly with step size. 

Let us consider terms up to in order to device a more accurate integrator. 



This equation contains the rate of change of acceleration and therefore is not very useful in this 
form. To simphfy these equations without loosing accuracy, let us specify j as 




(6) 



Vi{t + e) = Vi{t) + eai{t) + 0{e^) 
Xi{t + e) = Xi(t) + evi{t) + 0{e^) 



(7) 



Xi{t + e) = Xi{t) + evi{t) + -e^Giit) + 0{e^) 
Viit + e) = Vi{t) + eaiit) + ]^e^ji{t) + 0{e^) 



(8) 



m = 



a{t + e) - a{t) 



+ 0{e) 



(9) 



e 



This reduces the above equations to the following form 




(10) 
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e 


( Ax)}-m,s 


0.1 


5.4 X IQ- 


-4 


0.05 


7.3 X 10" 


-5 


0.02 


3.5 X 10- 


-5 


0.01 


1.2 X IQ- 


-5 


0.005 


5.3 X 10" 


-6 



Table 1: This table lists the values of e used in forward and backward integration of trajectories 
of a set of particles in an external potential and the corresponding error in recovering the initial 
conditions. Here we have used one grid length as the unit of length. 

This method can be used for velocity independent forces and is identical to the Leap-Frog method. 
The standard Leap-Frog method involves updating the velocity and position at an offset of half a 
step. To see the equivalence of the above equations and the Leap-Frog method let us consider the 
expressions for velocity at t + e/2 and t — e/2. 

v^{t + e/2) = Vi{t) + |ai(t) + ^e^ji{t) + 0{e^) 

Vi{t - e/2) = Vi{t) - ^ai{t) + ^e^ji{t) + 0{e^) (11) 

These two equations can be combined to give 

v^{t + e/2) = v^{t - e/2) + ea^{t) + 0{e^) (12) 
We can also use the expression for velocity at (t + e) to write 

Xi{t + e) = Xi{t) + €Vi{t + e/2) + ©(e^) (13) 
These two equations can now be combined to give the Leap-Frog method. 

Xi{t + e) = Xi{t) + evi{t + e/2) + O(e^) 

Vi{t + 3e/2) = Vi{t + e/2) + eai{t) + 0{e^) (14) 

This is called the Leap-Frog method as it updates velocities halfway between the step that is used 
to update the position. For velocity dependent forces these methods have to be modified into a 
predictor-corrector form. 

It is possible to improve the accuracy by using integrators that retain more terms in the Taylor 
series expansion. However most of these schemes require many evaluations of force for each step 
making them computationally uneconomical. 

2.3 Choosing the Time Step 

The time step for integrating the equation of motion can be chosen in many ways. In all of these 
methods one defines a suitable criterion [like energy conservation] and then obtains the largest time 
step that will satisfy the requirement. Following is a list of possible criteria : 
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• Monitoring validity of Irvine-Layzer equation. This is done by integrating the Irvine-Layzer 
equation [Irvine (1961); Layzer (1963); Dmitriev and Zel'dovich (1964)] 



— {2T + U) = C (15) 

a 

and monitoring the constant C. Variation in C with respect to total energy T + U can be 
used as an indicator of non conservation. [This is normally done with different forms of 
this equation so that only the kinetic or the potential energy term contributes to the integral. 
Estathiou et all d 1985^ ^1 In schemes where the force is computed at mesh points and then 
interpolated to particle positions the force at the particle position, f(x) ^ — Vxy?, i.e., the 
force does not equal the gradient of potential at a point. This leads to an apparent non con- 
serv ation of en ergy. If we include a correction term in the energy equation to account for this 
fact ( Efstathiou et all (Il985l) ) then it can be shown that it is possible to improve conservation 



of energy. This correction involves a direct sum over all pairs of particles and is therefore an 
almost impossible task to implement for a simulation with a large number of particles. Thus 
this is not a practical method of deciding the time step. 

• Convergence of final positions and velocities. If we evolve the trajectories of a set of particles 
with different time steps, then we expect the final positions and velocities to converge towards 
the correct value as we reduce the time step. This can be used to decide the time step. 

• Reproducibility of initial conditions. This is a more stringent version of the test outlined 
above. This method ensures that the solution we get is correct and errors are not building up 
systematically. If we run the particles forward and then back again we should in principle get 
back the initial positions. Although we ignore the decaying mode while integrating trajecto- 
ries forward in time, it tends to dominate if we evolve the system in the opposite direction. 
We can overcome this difficulty by not evolving the potential during this test. 

We find the last test to be the most stringent one and we use it to fix the time step for an arbitrary 
potential by the following construct. Consider an arbitrary potential ijj which has the maximum 
value of the magnitude of its gradient as Qmax = \^'^\max- We can write the equation of motion 
eqn.© for a small interval as 

e _ 24 

J^ffY ~ ■^S'maa; (16) 

Here we have used e to represent the largest displacement in time A0. We have specialised to 
r^o = 1 in writing the above equation. In which case, 9 = —2/y/a. For a given value of e we can 
fix the time step as 

A9 = 9^,P^^ (17) 



24g„ 

We can use the test mentioned above to fix the value of e. Table 1 lists a few values of e and 
{Ax)rms for a CDM simulation [boxsize=64/i^^Mpc] with 64^ particles in a 64^ box. The system 
was evolved from a = 0.02 to 0.25 and then back to a = 0.02. We use e = 0.005 for our simulations 
as we find the error in this case to be acceptable. 
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3 Calculating the Force 



An N-Body code solves the equation of motion and the Poisson equation in a self consistent manner. 
Therefore, moving the particles and computing the force for a given distribution of particles are the 
two most important components of an N-Body code. In this section we will discuss computation of 
force at mesh points for a given distribution of particles. For simplicity we shall assume that all the 
particles have the same mass. Generalisation to particles with different masses is not difficult but 
can strain the resources like RAM. [Setting up initial conditions for simulations that have particles 
with different masses is much simpler than setting up initial conditions for simulations that use 
particles with equal mass. See §4 for details.] 

In particle mesh codes the force at the particle position can be computed from the potential 
defined at mesh points in two ways. These are : 

• Energy Conserving Schemes : The force is computed at particle positions by using the gradi- 
ent of the interpolating function. 

f (x) = - Vx W{^, Xi)V'(xi)^ (18) 

The sum is over all the mesh points. An additional requirement is that the same interpolating 
function W should be used to compute the density field on the mesh from the distribution of 
particles. These are called energy conserving schemes as the force and gravitational potential 
are related to each other as f (x) = — Vx^(x). 

• Momentum Conserving Schemes : The force is computed on the mesh and then interpolated 
to the particle position. 

f (x) = - ^ W{x, (19) 

Momentum conserving schemes also require the same interpolating function to be used for 
computing density field on the mesh. It can be shown that in these schemes the force due to 
the ith particle on the jth particle is same as the force due to jth particle on the ith particle. 
[fij = — fjj.] This clearly is a necessary condition for momentum conservation. 

Use of mesh leads to anisotropy in force at scales comparable to the mesh size. A simple method 
for limiting these anisotropics is to choose a configuration that leads to a softer force at the mesh 
scale. Colli sionless evolution a lso requires the force to drop rapidly below the average inter-particle 



separation dMelott et al.l (119971) '). We choose to work with the momentum conserving schemes as 
these lead to a softer force at small scales in comparison with the energy conserving schemes. 
The algorithm for computing force at particle positions involves the following major steps. 

• Computing the Density Field : Mass of particles is assigned to mesh points by using an 
interpolating function. Density contrast is calculated from this "mass field" by using the 
definition 

M 
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Here M is the average mass per mesh point. Then fast Fourier transform technique is used to 
compute the Fourier components of density contrast. 

• Solving the Poisson Equation : The Poisson equation is solved in Fourier space. The FFT can 
be used to compute the gravitational potential on the mesh. 

• Computing the Force : The gravitational force equals the negative gradient of the gravitational 
potential, therefore the next step in algorithm is computing the gradient of potential. This 
can be done either in Fourier space or in real space by using some scheme for numerical 
differentiation. Force at the particle positions is computed by interpolation from the mesh 
points. 

In the following subsections we will discuss each of these steps in greater detail. 
3.1 Computing the Density Contrast 

Smoothing [interpolating] functions are used to assign mass to the lattice points for a given config- 
uration of particles. This is done by computing the following sum for every mesh point. 

M(xi)= J2 MparticleW{j^,^i) (21) 

all particles 

Here Xj are the co-ordinates of mesh points, x are the coordinates of particles and W is the interpo- 
lating function. As we have assigned same mass to all the particles we can scale the mass at mesh 
cites with the particle mass. 

Density contrast is computed from the mass defined at mesh points by using eqn. (l20t . FFT is 
then used to compute its Fourier transform. 

The smoothing function W should satisfy certain algebraic, computational and physical require- 
ments. 

• The sum of weights given to lattice points should equal unity. The interpolating function 
should be continuous and preferably differentiable as well. 

• The smoothing function should be peaked at small separations and should decay rapidly at 
large separations. The smoothing function should have the smallest possible base, i.e., the 
number of lattice separations over which it is non-zero should be small. The number of 
arithmetic operations required for interpolation is proportional to the cube of the number of 
lattice cells over which the interpolating function is non-zero. 

• The interpolating function should be isotropic. 

It is clear that we can not satisfy all these criteria fully and at the same time and therefore we have 
to consider a compromise solution. The first simplification that is made is to break up the three 
dimensional interpolating function into a product of three one dimensional interpolating functions. 
This reduces the complexity of the problem but also implies that the interpolating function cannot 
be manifestly isotropic. 
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These considerations, and the above mentioned simplifying assumption then leads to the fol- 
lowing set of interpolating functions. 

• Nearest Grid Point [NGP] : This is the simplest interpolating function. In this assignment 
scheme all the mass is assigned to the nearest grid point. 

where L is the grid spacing. This clearly satisfies the algebraic condition that the sum of 
weights should equal unity. However, this function is neither continuous nor differentiable. 
The leading term in the error in force of a point mass computed using this mass assignment 
scheme has a dipole like behaviour. 

• Cloud In Cell [CIC] : This the most commonly used interpolating function. Linear weights 
are used and mass is assigned to the two nearest grid points. The function is continuous but 
not differentiable. 

, . { l-\x-Xi\/L (for |x - Xj| < L) 

W{X,Xi) = < ^ ' (23) 

^ ' ^ I (for |x -Xi\> L) 

• Triangular Shaped Cloud [TSC] : The third function we wish to discuss is a quadratic spline 
that is continuous and also has a continuous first derivative. In this case the mass is distributed 
over three lattice points in each direction. 

[ I~{^y (forAx<L/2) 
W{x, Xi) = l I _ I ^ 1^ ' (for L/2 <Ax< 3L/2) (24) 
[ (for3L/2<Ax) 

Here Ax = |x — Xi| is the magnitude of separation between the particle and the mesh point. 
This function satisfies the requirements from the physical point of view but it is generally 
considered expensive from the computational point of view. 

It is obvious that these smoothing functions become more and more computationally intensive as we 
go for a larger base. All of the functions described above satisfy the essential algebraic requirements, 
so that the final choice is a compromise between the accuracy of the force field obtained using these 
and the computational requirements. However, we do not go beyond TSC as the base becomes very 
large and the transverse error in force becomes unacceptably large. 

Figure 2 compares these smoothing functions in a realistic situation. We have plotted the average 
weight assigned to a point that is at a distance r from the mesh point of interest. The curve in each of 
the panels shows the average weight where the average is computed over all directions [6, (j)] while 
keeping r fixed. The vertical bars show the rms dispersion about this average to depict the level of 
anisotropics in each of these functions. It is clear that TSC is the most isotropic of these functions 
and CIC is a "reasonable" compromise between isotropy and computational requirements. We will 
use the TSC smoothing function in the N-Body code for all applications. 
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Figure 2: This figure shows the weight assigned to a point that is at a distance r from the mesh 
point of interest. The curve in each of the panels shows the average weight where the average is 
computed over all directions [0,0], keeping r fixed. The vertical bars show the rms dispersion about 
this average to depict the level of anisotropics in these functions. 



12 



3.2 Solving the Poisson Equation 

The equation we wish to solve is 

V^^ = - (25) 
a 

This equation is to be solved with periodic boundary conditions. In particle mesh codes the equation 
is solved in the Fourier space. To solve this equation with periodic boundary conditions we must 
ensure that only the harmonics of the fundamental wave number kf = 2tt/NL contribute to the 
function. {NL is the size of the simulation box.] This type of sampling of k space is also required 
by the fast Fourier transforms. Therefore as long as we use FFT for computing Fourier transforms 
the correct boundary conditions are automatically incorporated. 

Two types of Green's functions are commonly used for solving the Poisson equation in Fourier 
space. These are 

• Poor Man 's Poisson Solver : This method uses the continuum Green's function. 

a 

G{k) = -^ (26) 

The main advantage of this is that it is the true Green's function. Also, it is isotropic by 
construction. 

• Periodic Green's Function : An alternative is to substitute Fourier transform of a discrete 
reaUsation of the Laplacian. A Green's function derived in this maimer tends to increase 
power at small scales and leads to larger anisotropics in the force field. A commonly used 
Green's function of this class is 

^^^^ ~ ~4 [sin^ k^L/2 + sin^ kyL/2 + sin^ k^L/2] ^^'^^ 

Here L is the grid spacing and equals unity in the units that we are using here. This Green's 
function is invariant under the transformation k^ ^ k^ + 27r. All FFT routines use such a 
transformation to rearrange wave modes for ease of manipulation [See §3]. Therefore this 
Green's function is easier to implement as compared to the continuum version as one does 
not have to worry about the rearrangement of wave modes. The main disadvantage of this 
Green's function is that it is anisotropic at small scales. 

We have plotted the two Green's functions in figure 3. Like figure 2 the curve here shows the 
average over angular coordinates for any value ofk = |k|. The vertical bars show the rms dispersion 
about this average and are an indicator of anisotropy. It is clear from this figure that the periodic 
Green's function enhances power at small scales, beyond the true value, and also introduces some 
anisotropy in the potential. 
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Figure 3: This figure shows the two Green's functions described in the text. Upper panel shows the 
continuum Green's function and the lower panel shows the Green's function derived from discrete 
representation of the Laplacian. Curves in each panel show the average of G{k) over angular coor- 
dinates for any value of A; = |k|. The vertical bars show the rms dispersion about this average and 
are an indicator of anisotropy. 
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CT 

Figure 4: This figure shows the rms error in solving Poisson equation as a function of the parameter 
a defined in the text. The dashed line shows the error for the periodic kernel whereas the thick line 
shows the error for the poor man's Poisson solver. The comparison was carried out on the grid so 
as to avoid mixing errors in solving the equation with errors introduced by interpolation etc. 
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We compared these two Green's functions by solving the Poisson equation for the density con- 
trast 

3 \ 




4 - ^ ) 



(28) 



2(72 

The solution to the Poisson equation for this source can be obtained analytically and the potential is 



tp = — exp 
a 



"2a2 



+ const. (29) 



We solved the Poisson equation numerically using the two Green's functions defined above on 
a 64^ mesh. We computed the rms deviation from the true solution in each case for a large range of 

values for a. The result of this analysis is presented in figure 4 that shows the error as a function 
of a. This curve clearly shows that the periodic kernel introduces more error than the poor man's 
Poisson solver. 

3.3 Computing the Gradient of Gravitational Potential 

Numerical differentiation is used to compute the force [— V^] on the mesh. The gravitational 
potential V' is also defined on the mesh. Interpolation is used to compute the force at the particle 
position. 

Here we will discuss three most commonly used methods for computing the gradient of potential 
in N-Body codes. 

• The most accurate and fastest method for computing derivatives uses fast Fourier transforms. 
The Fourier components of force are computed directly from the Fourier components of grav- 
itational potential. Only possible source of errors here is the FFT routine itself. The FFT 
routine used here gives an error of about 10^^% for Fourier transforming a sine wave in a 64^ 
box. [This error was estimated by comparing numerical values of Fourier components with 
their expected value.] 

The expression for Fourier transform of a derivative is 

f'k = -ikfk (30) 

• The simplest method for computing the derivative in real space is the mid point rule. 

+ Oih^) (31) 



df _ f{x + h)-f{x-h) , ^,^2^ 



dx 2h 
The equivalent expression for computing the derivative on the mesh is 

dl _ /(n + l)-/(n-l) 

dn ~ 2 ^ 

For an individual sine wave, the fractional error varies as fc^ for small k. Here, as for every 
operation in cosmological simulations, periodic boundary conditions are used to compute the 
derivative near the edges. 
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Figure 5: This figure sliows rms error in numerical differentiation. We differentiated sine waves 
with frequencies kx, ky and kz where ki = 27r/xj. We then compared the numerical result with 
its analytical counterpart to obtain an estimate of error. We have plotted the error as a function of 
X = 2t^ /{k1 + ky + k^Y^'^. Dashed line shows error for the five point formula and the thick line 
shows error for the mid point rule. Horizontal dot-dashed lines show the 10% and the 5% error 
level. This error analysis was carried out in a 64^ box. 
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Method 


5% error {k, I) 


10% error {k, I) 


Mid-Point 


0.7, 9.0 


1.0, 6.3 


Five Point 


1.4, 4.5 


1.7, 3.7 



Table 2: This table lists the scale at which error in differentiation crosses the 5% and 10% level for 
the mid-point and the five point methods. A comparison of these values suggests that the five point 
method is more robust and less error prone of the two. 



• It is possible to construct a more accurate formula by using the method of undetermined 
coefficients [See for example Antia (1991)]. A five point formula for differentiation that 
gives fractional errors proportional to A;^ for small k is 

dl ^ fix 2/0 - fix + 2h) + S{J\x + h) /{:,: h)) 
dx 12h 

+ Oih^) (33) 

This translates into the following expression for a uniformly spaced grid : 

df fin - 2) - fin + 2) + 8(/(n + 1) - /(n - 1)) 



dn 12 



(34) 



To compare these methods we studied the variation of relative error for sine waves as a function of 
frequency. We computed the numerical derivatives of a set of sine waves in a 64'^ box and compared 
these with the analytical results. We have plotted the root mean square error as a function of the 
inverse of wavenumber x = 2iT/kva figure 5. Dashed Une shows error for the five point formula 
and the thick line corresponds to the mid point rule. It is clear from these curves that the five point 
formula is far more accurate than the mid point rule. 

We have listed in table 2 the wavenumber for at which the average error crosses some thresholds 
for these two formulae. 



3.4 Tests of Force Calculation 

In the preceding subsections we have studied the steps involved in computing the force field for 
a given distribution of particles. For each step we discussed alternate methods of computing the 
required quantity and also discussed their relative merits. We shall now discuss the errors in com- 
puting the force field for one particle. We will use many combinations of individual components 
Uke the interpolating function, Green's function and the differentiator and compare the errors for 
each combination. 

To compute the error we placed a particle at a random position in a given lattice cell and com- 
puted force due to this at a fixed distance r in many directions to obtain the average force. We 
also computed the rms dispersion about the average. This process was repeated for one hundred 
positions of the particle in the lattice cell. Errors in this force were computed by comparing with 
the true 1/r^ force expected in such a situation. The error was split into two parts : transverse and 
radial. The results of this study may be summarised as follows : 
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Figure 6: This figure sliows tlie error in force of one particle. The top panel shows the root mean 
square dispersion in the radial force. The lower panel shows corresponding curves for the transverse 
component of force. The thick line shows the error for TSC interpolating function, the dashed line 
for CIC and the dot-dashed Une corresponds to NGP. In these panels the potential was computed 
using the poor man's Poisson solver and the gradient of the potential was computed in the Fourier 
space. 
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Figure 6: Continued. In these panels the Poisson equation was solved using the periodic Green's 
function and the mid point formula was used for computing the force. Notice that the CIC and TSC 
curves converge in this case indicating that the error is dominated by some other component. 
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The average force deviates very little from the true force for distances larger than one grid 
length [r > L] for almost all combinations of components. 



The root mean square dispersion about the average can be large for some combinations of 
components. 

Smallest rms dispersion is obtained for the combination of TSC, poor man's Poisson solver 
and numerical differentiation in Fourier space. 

• The combination of CIC, periodic kernel and mid point differentiation also gives small errors. 
Surprisingly the errors in CIC and TSC [other components remaining same] converge at scales 
larger than the grid size. This clearly indicates that the error is dominated by some other 
component, probably the mid-point differentiation. [Convergence disappears when five point 
formula or differentiation in Fourier space is used.] 

To elucidate some of the points made above we have plotted the radial and transverse error for the 
two combinations of components mentioned in the last two points. Figure 6 shows that the error is 
minimum for TSC amongst the three interpolating functions in all the cases. Other points mentioned 
in the above discussion are also brought out clearly in these panels. We will use the combination 
of TSC interpolating function and the poor man's Poisson solver. We will compute the gradient of 
potential in the Fourier space. 



4 Fast Fourier Transforms 



In previous sections we discussed methods used for moving the particles and the steps involved in 
computing the force in N-Body simulations. We mentioned that fast Fourier transforms are used for 
computing the Fourier transforms and play a vital role in reducing the time taken in computing the 
force. In this section we will start with the usual expression for the Fourier transform and reduce 
it to the form in which it is processed by the FFT For details o f the FFT technique s we refer the 
reader to any standard book on numerical analysis. ( Antial l 1991 >: Press et al. ( 198^) 1 



For simplicity we will work in one dimension in this section. However we will write the d di- 
mensional generalisation of the main result. Consider a function f{x) that has the Fourier transform 
g{k). Then we may write 



fix) 



■ikx 



(35) 



It is not possible to use computers to integrate arbitrary functions over an infinite range in finite 
time. The integral is thus truncated at some finite kmax- This truncation is possible only for those 
functions for which the integral over \k\ > kmax is sufficiently small. Thus the integral to be 
evaluated becomes 



fix) 



^9ik)e-^'^^ 



(36) 
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This integral is evaluated in form of a summation over the integrand. 



where the function is sampled at intervals of AA; in this interval. For a constant Afe, this sampling 
will divide the range —kmax to kmax in N intervals where N = {2kmax/^k) — 1. The fast Fourier 
transforms require the number of intervals to equal an integer power of 2, which clearly is not 
possible for the above range if we are to sample the A; = mode as well. It is possible to circumvent 
this problem by making the range of integration asymmetric. We add one extra interval after kmax 
and thus make the number of intervals equal to an even number. We can further choose kmax and 
Ak to ensure that A?^ is an integer power of two. Thus the sum may be written as 

+Ak 

kmax 

1 N/2 

= ]^ E 5(Afc m)e-^^'=-- (38) 

m=-N/2+l 

The boxsize NL = 2-K/Ak defines the scale over which the function f{x) is defined apart from its 
periodic copies. Here N = {2kmax + Ak)/ Ak = 2^ where j is a positive integer. The function 
f{x) is defined on a regular mesh with L as the spacing between adjacent mesh points. This allows 
us to rewrite the previous equation as 

N/2 

f{nL) = — Yl 5(AA; m)e-^2Wiv (39) 

m=-N/2+l 

To rewrite this in the form used by the FFT algorithm we split this sum into two parts : 

f{nL) = ^ 5(AA;m)e-2'^-"/^ 

m=-N/2+l 



+ ^^^^ m)e-^2Wiv (40) 



N/2 



m=0 



Changing the summation index in the first sum hy N [m ^ m + N] and dropping a 2-jt factor in 
the phase, we obtain 



N-l 



f{nL) = ^ E 5(Afc(m-iV))e-2™„/iV 



m=N/2+\ 

N/2 



NL ^ 

m=u 
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These terms can be put together in the following form 

, N-l 

finL) = ^ E Gi^k m)e-^2WiV ^42) 

m=0 

where G{Ak m) = g{Ak m) for m < N/2 and G(AA; m) = g{Ak{m-N)) for m > N/2. From 
here it is clear that the input array for FFT should have the negative wave numbers after the positive 
wavenumbers. [This can also be interpreted as "periodicity" in k space.] The generalisation of the 
above expression to d dimensions is 

f{nL) = — ^ J: G(AA;m)e-2'^'"-/^ (43) 

^ ' m=0 

Here G{m) is the Fourier transform of the function f{nL) with every component rrii > N/2 
replaced by rrii — N. 

If an FFT routine uses a different normalisation then the input array must be rescaled to get the 
correct ampUtude for the function /(nL). If the FFT routine uses 

1 

/(nL) = j^Yl G'(AA;m)e-^2'^™-"/^ (44) 

^ m=0 

then we must scale the input function as G^^'^ = G* {l/NLy. This scaUng is not important if we 
are computing both the forward and the inverse transform as the overall normaUsation must be same 
for all FFT routines. This scaling becomes relevant only in those cases where we are transforming 
a given quantity only once, as in the generation of initial potential. 

The above equations estabUsh the form in which an array must be passed to the FFT routines in 
order to compute the Fourier transform. 



5 Setting up Initial Conditions 
5.1 The Initial Density and Velocity Field 

N-Body simulations are generally started from fairly homogeneous initial conditions, i.e. the density 
contrast is much smaller than unity at all scales that can be studied using the simulation. In this 
regime we can use linear theory to compute all quantities of interest. 

In linear theory, the evolution of density contrast can be described as a combination of a growing 
and a decaying mode. We can write the solution for perturbations in an Einstein-deSitter Universe 
as 

(5(a) = cia + C2a~^/^ (45) 

We can evaluate ci and C2 by using the initial conditions. These conditions can be written in terms 
of the initial velocity field and the initial potential. 

iT) = -^-"^ (46) 

\oaJ in 
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Using these we can express the linear solution in terms of the initial potential and the initial velocity 
field as 



5(a) 



3^2 , 2^ 

5 5 



2 

«+5 



a"3/2 (47) 



This equation suggests that for a system in growing mode, Uj„ = —Vipin- This equation gives us 
the density contrast at each point in terms of the initial gravitational potential. [We will discuss 
the problem of generating the initial potential in the next subsection.] Given an initial potential the 
above equation suggests two schemes for generating the initial density field. These are : 

• The particle are distributed uniformly and their masses are proportional to the local value 
of density. We can either start with zero velocities, in which case we have to increase the 
amplitude of V'm by a factor 5/3 to account for the presence of decaying mode. Alternatively 
we can choose to put the system in the growing mode and assign velocity Ui„ = —Vijjin to 
each particle. One major drawback of this method is that an extra array containing masses of 
particles needs to be stored and this can be a problem for large simulations. 

• Starting with a uniform distribution the particles are displaced by a small amount, say much 
less than the inter-particle separation, using velocity Uj„ = —Vipin- The resulting distribution 
of particles will represent the required density field if the initial distribution did not have any 
inhomogeneities. We can also retain the initial velocity field. 

Thus the main requirement from initial positions of particles before we generate the required per- 
turbation is that the distribution of particles sample the potential "uniformly." Any inhomogeneities 
present in the initial distribution will combine with the density perturbations that are generated by 
displacing particles and will modify the initial power spectrum. 

One obvious solution that ensures a uniform distribution is placing particles on the grid. This 
generates a distribution that is uniform but not random. An additional problem with this distribution 
is the extreme regularity which leads to shot noise at Nyquist frequency for potentials with large 
coherence length when the number of particles is smaller than the number of cells in the mesh. [See 
figure 7] 

Another "obvious" solution is to put particles at random inside the simulation box. This distri- 
bution is uniform but it has \/]V fluctuations which result in spurious clustering. The fluctuations in 
number of particles in a cell of size R decrease as 1/ R^^^ with scale. [Comparing with power law 
spectra we note that this is an n = spectrum.] 

Apart from these simple minded solutions there exist [at least] two other distributions of parti- 
cles which may be used for the initial positions. 

The Glass initial conditions are obtained by running a random distribution of particles through 
N-Body with a repulsive force. In this case we start with small perturbations and the repulsive force 
leads to a slow decrease in the amplitude of perturbations. As the perturbations get progressively 
smaller it worthwhile to estimate the evolution of inhomogeneities using linear theory. The equation 
for evolution of perturbation s in this case is same as the standard equation for linear growth of 
perturbations ( Peebles ( 1980l) ') but with a source term that has an opposite sign. This leads to the 
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following oscillatory solution with a decaying amplitude. 



5ia) 



-1/4 



Q COS 



23, \ ^ . / V23, ^ 
ma + p sin — - — in a 



(48) 



Where a and (3 are constants to be fixed by initial conditions. The process of generating initial 
conditions can take a long time for a large number of particles. However this has to be done only 
once and the data can be stored in a file for repeated use. 

The other method is a variant of the grid and Poisson initial positions mentioned above. Here 
particles are placed in lattice cells but at a random point within the cell. This removes the regularity 
of grid without sacrificing uniformity. The number fluctuation falls faster than ^/N and have a 
smaller amplitude even at the smallest scale. The amplitude of fluctuations can be controlled by 
reducing the amplitude of displacement about the mesh point. 

After the initial density field has been generated the initial velocities for N-Body can be set as 
Uin = — VV'm if we want the system to be in the growing mode. However, generation of density 
field from the initial potential, as outlined above involves many numerical operations that modify 
the potential at scales comparable to the mesh size. A better method, for consistency of density, 
velocity and potential, is to recompute t he potential after gener ating the density field and set the 
velocities with the recomputed potential dEfstathiou et al.l (119851) '). This is particularly important in 
case of models with a lot of small scale power. If we use the initial potential for fixing velocities 
then velocities have a stronger small scale component than would be expected from the density field 
that has been generated, leading to inconsistency in the input configuration for N-Body. Therefore 
we use Uin = — VV', where V^V = S{a)/a is the potential generated from displaced distribution of 
particles. 



5.2 Initial Potential 

In this section we shall outline the method used to generate the initial potential. The same method 
can also be used to compute the density field directly if the initial density field is to be generated by 
placing particles with different masses on the mesh. 

In most models of structure formation the initial density field is assumed to be a Gaussian 
random field. Linear evolution does not modify the statistics of density fields except for evolving 
the amplitude. As the potential and density contrast aie related through a linear equation, it follows 
that the gravitational potential is also a Gaussian random field at early epochs. 

A Gaussian random field is completely described in terms of its power spectrum. The Fourier 
components of a Gaussian random field [both the real and the imaginary part] are random numbers 
with a normal distribution with variance proportional to the power spectrum of the random field. 
The proportionality constant depends on the Fourier transform convention. To fix this constant we 
consider the probability functional for a function g. 



P[g]=Bexp 



d^k |g(fc)p 
(2vr)3 Pg{k) 



(49) 
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Here S is a normalisation constant. In the Fourier transform convention used in §3 this equation 
can be written as 



P[g]=Bexp 



1 1 



(50) 



Thus the variance of the real and imaginary components for each Fourier mode of a Gaussian ran- 
dom field should equal 



(51) 



The factor of 2 takes into account the fact that both the real and the imaginary components share the 
variance. Thus the function g may be written as 



g{k) = {ak +ihk) 



.1/2 



(52) 



where and 5^. are Gaussian random numbers with zero mean and unit variance. 

To generate the gravitational potential we substitute the gravitational potential ipf^ in place of 
g{k) and the power spectrum with the power spectrum for the gravitational potential. The power 
spectrum for the gravitational potential is P^{k, a) = Ps{k, a) / {a?k^) = P^^^{a = \,k)/k^. Here 



can write 



5 yu — l,k) h the linearly extrapolated power spectrum of density fluctuations. With this, we 



pii-^(a =l,k)N^L^ 



2k* 



1/2 



(53) 



Here we have taken a (to) = 1- To specialise to the Fourier transform convention used in the FFT 
routine we use, we use / = iV^/^. reqn. (l44H This implies 



FFT 



Ar3/2^3 



(ofc +ibk) 



Ps{k) 



.1/2 



2L^k* 



(54) 



Here k = 2TTm/{NL). If Ps{k) is of the form Ak"-f{k) where / is a dimensionless function [It 
may be the transfer function or a cutoff imposed by hand.] and A is the amplitude. In the units we 
are using here L = I, therefore 



where all lengths are written in units of L. 



A /2iTm 



n-4 



2 \ N J 



/27rm\ 



1/2 



(55) 



5.3 Testing Initial Conditions 

In this section we shall test the initial conditions generated by the method outUned above. We will 
compare the power spectrum of density perturbations generated by this method with the theoretical 
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L=27T/k 

Figure 7: This figure sliows the theoretical power spectrum and the induced power spectra, linearly 
extrapolated to a = 1. The theoretical power spectrum is shown as a thick Une and the power 
spectrum generated by displacing the particles is shown as filled squares [for simulations with 128^ 
particles] and as empty circles [for simulations with 64^ particles]. The upper panel here corre- 
sponds to the power law model with n = — 2 and lower panel to the model with n = 0. 
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L=27T/k 



Figure 7: Continued. These panels show the initial power spectra for models with Gaussian power 
spectrum. The upper panel corresponds to a Gaussian power spectrum with peak at Lq = 8L and 
the lower panel corresponds to a spectrum with Lq = 24L. The secondary peaks at small scales 
occur at the harmonics of the wavenumber corresponding to the peak of the initial Gaussian. 
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power spectrum for a wide range of models. This comparison was carried out using simulations 
done with 128'^ particles in a 128'^ box. To study the effect of a smaller number of particles we also 
studied the power spectrum for a simulation with 64^ particles in a 128^ box. We use the following 
models for this purpose : 

• Power law models with n = — 2 and n = 0. 

• A Gaussian power spectrum with spread ak = 2it/NL peaked at /cq = 2it/Lq. \_P{k) oc 
exp{— (/c — A;o)^/2(T^}] This model was studied for Lq = 8 and 24. 

Models outlined in the first item have power at all scales whereas for the Gaussian power spectrum 
the power is concentrated in a narrow range of scales. This tests the accuracy of initial condition 
generation for the two extreme classes of models. 

We have plotted the theoretical power spectrum and the induced power spectra, linearly extrap- 
olated to a = 1, in figure 7 for the models listed above. The theoretical power spectrum is shown as 
a thick line and the power spectrum generated by displacing the particles is shown as filled squares 
for simulations with 128'^ particles and as empty circles for simulations with 64^ particles. 

The simulated power spectrum agrees with the theoretical power spectrum for scales larger than 
about 5-L for power law models. In this range the power spectrum in simulations with different 
number of particles also agrees quite well. 

In case of Gaussian power spectrum the simulated power spectra also contain some power at 
harmonics of the main peak, i.e., at scales Lq/2, Lq/3, etc. The Gaussian at scale Lq is reproduced 
quite well. Simulations with 64"^ particles contain spurious noise around the Nyquist scales. Apart 
from this difference power spectra in two cases agree quite well. 

Simulations with 64^ particles generally have more noise at smaller scales. For hierarchical 
models this is essentially shot noise. However for models with large scale coherence [as in case 
of Gaussian power spectrum with Lq = 24L or n = — 2 power law model] this can lead to large 
spurious noise at the Nyquist scale. This can be understood in terms of the large coherence length 
in the displacement field for such models, which leads to a sequence of alternately filled and empty 
cells. This contributes to the excess noise at the Nyquist scale. 

The theoretical power spectrum is reproduced correctly in a large range of scales for a wide va- 
riety of models. Power spectrum at scales smaller than about 5L does not match with the theoretical 
power spectrum and hence results in this regime should not be used at early times. This constraint 
applies to all types of N-Body codes as the same method is used to generate the initial conditions 
for all types of N-Body codes. This constraint does not apply at late times as the power at sm all 
scales at late times is dictated by initial power at larger scales. ( Ba gla and PadmanabhanI ( 1997 )) 



6 Tests of N-Body Code 

In this section we will test performance of the N-Body code as a whole. We will compare the 
evolution of gravitational clustering in N-Body simulations with analytical results. The number of 
analytical results available in the non-linear regime is very limited, and hence the number of tests is 
also very small. 
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k/kf 






1 


0.021 


0.034 


2 


0.048 


0.073 


4 


0.104 


0.138 



Table 3: This table lists the errors in positions and velocities for one dimensional collapse. The 
errors are computed by comparing with Zel'dovich approximation. The errors have been listed 
for three different initial perturbations, the first column contains the wavenumber[in units of the 
fundamental wavenumber kf = 27r/N] of the sine wave perturbation. The errors were computed at 
the time of first shell crossing in the system. The N-Body simulation used for this was carried out 
using 64^ particles. 



As a test of evolution of positions and velocities we can compare the motion for a one dimen- 
sional collapse before shell crossing with Zel'dovich approximation (Zel'dovich ( 1970)). This test 
was suggested by Efstathiou et al.(1985) and they compared the positions and velocities for such 
a collapse with the corresponding quantities in Zel'dovich approximation. They chos e a sine wave 
along one of the coordinate axes as the initial perturbation in potential. Recently (^ Melott et al. 
(12221)) this test has been generalised so that the plane wave is not along any coordinate axis. In 
such a case the particles can approach each other with a small nonzero impact parameter. Monitor- 
ing velocities along directions perpendicular to the perturbation provides a good test of collisionless 
evolution. Collisionless evolution of a one dimensional perturbation can not induce velocities along 
directions orthogonal to the perturbation. The authors find that Particle mesh codes are the only 
ones that pass this test. P^M and tree codes can also provide collisionless evolution if the softening 
length is of the same order as the inter particle separation. 

We carried out this test with one dimensional perturbation along the x axis. The results are listed 
in table 3. Here we have tabulated the root mean square error in displacement and velocity. These 
are defined here as 



Xi 



1/2 



1/2 



(56) 



The table lists errors for three wavenumbers. In each case we have given the error at the epoch of 
first shell crossing in the system. The error increases with the wavenumber of the sine wave used for 
initial potential as we sample a given wave with a smaller number of particles and the discretisation 
effects become more important. 

The above test checks the working of an N-Body code for a very special case. In a more general 
situation it is difficult to compare the output of numerical simulation at the level of positions of 
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Figure 8: This figure siiows h{a,x) as a function of ^{a,x) for two power law models [n = 
and —1] and the CDM model. The thick line depicts the relation between these quantities in linear 
theory and the results from simulations are shown as points. All the points in the range ^ ^ 1 lie 
along the line with little dispersion. 



particles and only a statistical comparison is possible. We will use the averaged correlation function 
^ and the scaled pair velocity h (Peebles ( 1980)) for testing the accuracy of the N-Body code. 

In linear regime the averaged correlation function f evolves as a^, i.e., ^ (x, a) = {a^ /af)^{x, ai) 
In addition we know that £ ^ 1 in this regime. These can be used, along with the pair conservation 

I 1 J 1 D JT 

equation Peebles (1980; to show that /i = 2^/3 for £ ^ 1. This equation is valid for all models 
in linear theory. We have plotted /i(a, x) as a function of ^(a, x) in figure 8. We have plotted the 
line corresponding to the linear theory result and have shown results from N-Body simulations as 
points. We used simulations of power law models [n = and —1] and the standard CDM model 
for this plot. Points corresponding to different simulations have been shown with different symbols. 
All the points crowd around the line h = 2^/3 in the linear regime with little dispersion. The 
points deviate from this line for large f but the dispersion between different models remains small. 



This c o rresponds to an approximate "universa lity" in the relation between h and £. IHamilton et al. 
( 199lb :lNitvanan da and Padman abhan'(' 19941)1 



Self similar evolution of £ for power law models (iPeeblesI (119801) ') can be used to test correctness 
of the non-linear evolution of gravitational clustering. In figure 9 we have plotted ^ (a, x) as a func- 
tion of X /xni for the n = power law model. Here Xni is the scale where the linearly extrapolated 
^ is unity. This scale is proportional to a?/^ for the n = power law model. If the evolution is 
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Figure 9: This figure shows ^(o, x) as a function of x jx^i for the n = power spectrum. The points 
have been plotted for two epochs. The overlap between points clearly shows that the system is 
evolving in a self similar maimer. The dashed line shows the linear slope of the averaged correlation 
function. 
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Figure 10: Comparison of simulations at different scales. This figure shows ^^/^ for three CDM 
simulations done with boxsize of 64, 128 and 256/i~^Mpc. Matching of curves in the overlapping 
region shows that the numerical simulation is not introducing any artifacts. The dashed line shows 
the same function from a P^M simulation by Brieur et al.(1995). This simulation was done with the 
same initial power spectrum but a different realisation of the Gaussian random field was used. The 
similarity between ^ in the three PM simulations and the P^M simulation shows that the non-hnear 
evolution is being followed in the same manner in both the simulations. 

self-similar then the points from different epochs must he along a single curve. Figure 9 clearly 
shows that the evolution of averaged correlation function follows a self similar pattern over a large 

range of scales. 

The last test we consider here compares the averaged correlation function for simulations of the 
standard CDM model. We compare ^^/^ for three simulations carried out with boxsize 64, 128 and 
256/7."^Mpc. Fi gure 10 shows that the curves match in the overlapping region, implying that the 
numerical simulation is not introducing any scale in the non-linear evolution of density perturba- 
tions. We also compare these curves with ^1"^ obtained from a different simulation. This reference 
simulation was done by Brieur, Summers and Ostriker (1995) using a P^M code implemented on a 
GRAPE machine. This simulation was done with the same theoretical power spectrum but using a 
different realisation of the initial density field. The two types of simulations being compared have 
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virtually nothing in common as far as the implementation of the mathematical model is concerned. 
The similarity in these curves implies that both the codes aie evolving density perturbations in the 
same manner. 



7 Analysis of N-Body Data 

In this section we will outline the method used for computing the correlation function and the scaled 
pair velocity (Peebles (198Q,) ) from simulation data. 
The averaged pair velocity h is defined as 



/i(a, x) 



/ «(^-^)-(x.-x,) \ 

El 



(57) 



In this equation subscripts label particles and the averaging is over all pairs with separation equal to 
X, i.e., x"^ = (xj — Xj)^. This quantity is computed from simulation data by binning the pairs by pair 
separation x. The number of pairs in a given bin as well as the quantity to be averaged is summed 
for each pair in each bin. The ratio of these quantities for each bin gives the value of h{a, x). 

The correlation function was defined in §1.4 as the Fourier transform of the power spectrum. 
This is equivalent to the following definition 

aa,x) = {6{y)5{z)) ;x = |y-z| (58) 

where the average is over all pairs of points with separation x and over all x with the same magni- 
tude. This can be rewritten as 



1 (59) 



g(a,y) g(a,z) 
Q{a) g{a) 

Here we have used the fact that density contrast is a random field with zero mean. We can replace 
the densities by number densities and the product of number densities at points separated by a given 
distance by the number of pairs with that separation. With this the above equation reduces to 

C{a,x) = ^^-1 (60) 
n{a) 

where n(a, x) is the number of pairs with separation x and n is the number of pairs in a uniform 
distribution. Therefore the problem is again reduced to computing the number of pairs separated by 
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a given distance. This is done by dividing the range of x into small intervals and binning the pairs 
into these intervals. Similar operation is required for computing the scaled pair velocity, therefore 
the binning operation for computing these quantities can be combined conveniently. 
We can obtain an expression for the averaged correlation by averaging ^. 

(,{a,x) = j dyy i[a,y) = 

X 

I dynjy) 

= ^ 1 (61) 

/ dyn{y) 

where we have used the fact that the average number of pairs with a given separation is proportional 
to the square of the separation for a uniform distribution of particles. In the discrete realisation of 
this expression the integrals over number of pairs are replaced by summation over bins. Compu- 
tationally this quantity is less error prone than the correlation function as for large separations the 
number of pairs in a given bin in the simulation output approaches the number of pairs for a uni- 
form distribution. Subtracting two nearly equal numbers can give large error. However, in case of 
the averaged correlation function the excess number of pairs at small scales is carried over to large 
scales through summation over all smaller scales. And in general the ratio {J2 ^) I ^) larger at 
all scales in comparison with n/n. 



8 Discussion 

In the preceding sections we have described the basic mathematical model that is implemented in 
particle mesh cosmological N-Body codes. We have not discussed the P'^M c odes, tree codes, etc . 
in this review as only PM codes are known to ensure coUisionless evolution (iMelott et al.l dlQQTn i 
though they suffer from a very limited resolution. Other methods improve spatial resolution but do 
not ensure coUisionless evolution. However, most of the the machinery is common to these codes 
and many of the results can be carried over to these codes. 
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